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ABSTRACT 


The  dimensionless  form  of  the  conduction  equation  pre¬ 
sented  in  this  thesis  is  mathematically  linear  and  has  been 
applied,  in  the  past,  to  several  fixed  boundary  conditions. 
With  phase  change,  the  overall  problem  is  non-linear  in 
nature.  The  reason  for  this  is  that  one  boundary  is  moving 
continuously  and  its  position  is  not  known  a  priori. 

A  theoretical  and  experimental  investigation  is  pre¬ 
sented  for  freezing  or  melting  in  one  dimension  when  the 
boundary  temperature  varies  with  time.  An  explicit  solution 
and  an  alternative  quadrature  solution  is  given  for  the 
semi-infinite  domain  when  the  boundary  temperature  varies 
as  a  simple  power  of  time.  Experimental  results,  for  an 
ice-water  system,  are  presented  for  particular  values  of 
this  power  and  are  compared  with  predictions  from  the 
quadrature  solution.  The  simple  power  law  solution  has 
been  extended,  in  approximate  form,  to  predict  depth  of 
freezing  and  melting  during  one  cycle  of  a  sinusoidual 
surface  temperature.  Comparison  with  experimental  results 
shows  qualitative  agreement  and  furnishes  new  insight  into 
the  problem  of  periodic  freezing  and  melting. 

Experiments  in  two-dimensional  ice  formation  are  in¬ 
cluded.  Typical  ice  profiles  are  presented  for  the  case 
of  two  surfaces  at  different  uniform  temperatures  with  one 
surface  held  below  the  freezing  temperature. 
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NOMENCLATURE 


x 


9/  f,  f 

T  ,  t 
X 

(3,  b 
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e 
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S 

P,  Q, 

a,  A,  B,  C 

K 

C 

P 

L 

P 

/c 

r 

CD 

a 

i  i  i 

r 


o 

i 

in ,  c 


temperature 

reciprocal  temperature  constant  as  defined 
time 

distance  from  domain  edge 
interface  depth 
similarity  variable 
dummy  variable 

temperature  index  (power  exponent  of  time) 
Laplace  transformation  variable 
functions  of  t)  defined  in  text 
constants 

thermal  conductivity 
specific  heat 
latent  heat  of  fusion 
density 

thermal  diffusivity 
gamma  function 
angular  velocity 
radians 

first  and  second  derivatives  (with  respect  to  t)) 
Laplace  transform  with  respect  to  time. 

SUBSCRIPTS 
location  x  =  o 
initial  condition 
maximum 


r_T  •••  -2E  ;  1 


**•»  -a 


- 

'  . 


- 


. 


CHAPTER  I 


INTRODUCTION 


When  a  liquid  undergoes  an  unsteady  change  to  the  solid 
phase,  thermal  energy  in  the  form  of  latent  heat  of  fusion 
is  released  at  the  moving  liquid-solid  interface,  the 
history  of  which  is  unknown  a  priori.  This  latter  pecu- 
larity  gives  the  problem  a  mathematically  non-linear  char¬ 
acter  and  relatively  simple  analytic  solutions  are  dif¬ 
ficult  to  obtain. 

More  than  a  century  ago,  Neumann  (Carslaw  and  Jaeger 


)  found  an  exact  solution  relating  to  the  solidifica¬ 
tion  of  the  semi-infinite  region  of  liquid  initially  above 
or  at  the  stable  crystallization  (freezing)  temperature 
with  the  boundary  temperature  suddenly  decreased  to  and 
maintained  at  a  temperature  below  freezing.  Stefan 
(Ingersoll  3  )  later  published  the  solution  for  the  simpler 
situation  when  the  region  is  initially  at  the  freezing 
temperature.  The  literature  frequently  refers  to  the 
subject  as  the  Stefan  problem. 

With  increasing  interest  in  the  casting  of  metals. 


Lightfoot 


found  an  approximate  solution  for  the  case 


of  the  slab  of  finite  thickness,  using  the  Neumann  boundary 


*  Numbers  in  brackets  denote  references  given  on  page  69  . 
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conditions?  it  was  shown  that  the  velocity  of  the  inter¬ 


face  was  practically  constant  upon  approaching  the  centre 
of  the  slab„  As  a  comparison,  this  thesis  contains  the 
results  of  an  experiment  using  a  slab  insulated  on  one 
surface  and  with  a  time  dependent  temperature  condition 
on  the  other o  A  more  recent  extension  was  given  by 


Kreith  and  Romie 


who  considered  a  constant  fusion 


velocity  for  cylinders  and  spheres  ?  in  addition  they 
considered  the  semi-infinite  domain  and  the  equation  for 
interface  depth,  (3 ,  in  a  dimensionless  form  which  is  id¬ 
entical  to  the  equation  termed  the  Neumann  solution  in 


this  thesis 


In  1950,  Landau 


formulated  the  governing  equ¬ 


ations  in  general  form  and  used  a  finite  difference  pro¬ 
cedure  to  present  a  set  of  solutions  for  the  semi-infinite 
region?  this  work  has,  in  the  past,  been  used  as  a  stand¬ 
ard  for  comparison  purposes „  In  1959,  Murray  and  Landis 
7  extended  the  range  of  initial  and  boundary  conditions 
with  a  numerical  approach  permitting  a  continuous  deter¬ 
mination  of  the  fusion  front  travel  and  indicating  the 
internal  temperature  distribution  with  good  accuracy. 

Several  diversified  approaches  have  been  presented 
in  solving  for  specific  boundary  and  initial  conditions  * 


Goodman 


8 


introduced  the  heat  balance  integral  following 


the  momentum  integral  principle  of  von  Karman.  The  concept 
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makes  use  of  a  second  order  polynomial  to  describe  an 
assumed  temperature  distribution  in  the  solid  region?  if 
the  assumption  is  correct,  accurate  analytical  solutions 


■will  resulto  Poots 


9,10 


used  integral  methods  to  sat¬ 


isfy  the  boundary  conditions  for  one-  and  two-dimensional 
phase  change  problems .  The  method  may  be  useful  for 
theoretical  comparison  with  the  two-dimensional  experi¬ 
ments  given  in  this  thesis. 

The  mathematical  generalization  of  the  Neumann -Stefan 
problem,  assuming  the  boundary  temperature  as  an  arbitrary 


function  of  time,  was  originated  by  Portnov 


11 


in  1962 


The  initial  temperature  distribution,  an  arbitrary  func^ 
tion  of  distance,  is  included  in  the  Laplace-integral 
solution o  The  method  attacks  the  problem  directly  but 


Westphal 


12 


shows  that  its  application  to  Arctic  ice 


formation  is  extremely  difficult 0  The  resulting  power 
series  appears  to  be  suitable  for  numerical  computer 
programming  for  various  boundary  conditions. 

The  purpose  of  this  thesis  is  to  present  an  "inter¬ 
mediate"  analytical  solution  extending  the  basic  work 
of  Neumann  and  yet  not  having  the  complexity  of  the 
Portnov  generalization.  The  solution  mathematically  de¬ 
scribes  the  solid  region  temperature  and  interface  lo¬ 
cation  in  dimensionless  form  employing  a  simple,  time- 
dependent  boundary  condition.  Only  the  case  of  the  liquid 
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initially  at  the  freezing  temperature  is  considered,.  The 
solution  is 

(1)  in  explicit  form  convenient  for  application  with 
standard  tables  and, 

(2)  in  quadrature  form  for  computer  calculation  of  inter¬ 
mediate  points  not  obtainable  by  the  explicit  solution. 

A  numerical  comparison  of  the  explicit  and  quadrature 
solution  is  included  as  Appendix  B.  An  application  of 
the  Laplace-integral  method  is  presented  as  Appendix  C 
for  comparison  with  the  explicit  solution  using  a  linear 
temporal  variation  in  boundary  temperature.  In  addition, 
the  case  of  the  sinusoidual  boundary  temperature  has  been 
considered?  an  approximate  solution  is  obtained  for  temp¬ 
erature  and  interface  depth  using  superposition  of  explicit 
solutions « 

An  explicit  solution  is  also  obtainable  for  the  case 
of  the  liquid  region  initially  at  a  uniform  temperature  above 
freezing  but  is  not  presented  in  this  thesis.  The  case 
of  an  initially  non-uniform  temperature  in  the  liquid  has 
not  been  investigated. 

Experimental  results  show  good  agreement  with  theory 
using  an  ice-water  system  and  simple  time-dependent  bound¬ 
ary  conditions.  Qualitative  agreement  is  shown  for  the 
sinusoidual  boundary  condition.  As  a  natural  extension 
of  the  unidimensional  problem,  experimental  results  for 
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the  two-dimensional  steady  state  case  have  been  included 
For  completeness,  a  brief  experimental  analysis  of  super 
cooling  is  presented* 
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CHAPTER  II 


ONE -DIMENSIONAL  THEORY 


2  o 1  FORMULATION  OF  THE  PHASE  CHANGE  PROBLEM 

Consider  the  semi -infinite  domain  (x  >  0)  shown  in 
figure  2.1.  Let  the  initial  temperature  of  the  liquid 
domain,  in  this  case  water,  be  at  the  melting  temperature 
which  shall  be  considered  the  temperature  datum.  The 
melting  temperature  is  used  as  the  fixed  reference  since 
the  " freezing”  temperature  is  variable  with  the  effect 
of  supercooling.  When  the  surface  temperature  9Q(t)  is 
lowered  then,  ideally,  ice  will  form  as  shown.  If  the 
ice  properties  are  assumed  to  be  constant,  the  differen¬ 
tial  equation  governing  heat  conduction  within  the  ice 
region  is  given  by 

c)^9  _  1.  <39  _  n 

9x2  "  *  St 

Mathematically,  then,  the  problem  consists  of  determining 
the  solution  to  equation  (1)  subject  to  the  initial  con¬ 
dition  9^ (x)  ~  0  and  boundary  conditions 

x  =  0  ;  9  =  9Q(t) 

x  =  b  (t)  ;  9  =  0 

where  b(t)  is  the  ice  depth.  Although  the  time  dependent 
boundary  condition  at  the  origin  may  take  many  different 
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forms  in  practice,  the  surface  temperature  will  be  re¬ 
stricted  to  the  form  0Q(t)  =  Ctn  where  C  and  n  are  con¬ 
stants.  Using  this  surface  temperature  variation,  a 
closed  solution  for  ice  depth  is  developed. 

Taking  the  similarity  forms  for  the  ice  region  as 


0(tj)  =  Ctn  f  (tj)  (3) 

n 

where  tj  =  x/2  ’V  tc  t , 

then,  equation  (1)  may  be  transformed  into  the  linear 
ordinary  differential  equation 


f '  '  +  2  T|f '  -  4nf  =  0 

with  transformed  boundary  conditions 

f (0)  -  1,  f(p)  -  0 

where  (3  =  b(t)/2  V  *  t. 


(4) 

(5) 


which  is  the  dimensionless  ice  depth. 


. 


.  V* 


' 
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2.2  SOLUTION  FOR  TEMPERATURE 


2.2-1  Explicit  Solution 

If  ip  is  a  solution  of  (4)  , 
may  be  written 

f  (t})  =  A  ip  +  B  ip 


then  the  complete  solution 
2 


or 


f  (n) 


ip 


A  +  BP  (r)) 

1 


(6) 


A  known  solution  for  ip,  using  the  surface  temperature 
©o(t)  =  Ctn,  is 

ip( rj)  =  22n  T(n+1)  i2nerfcr)  .  (7) 


Substituting  the  known  solution  (7)  into  equation  (6) ,  and 
satisfying  the  boundary  conditions  (5) ,  yields  the  tempera¬ 
ture  in  the  frozen  region  as 


where 


f  (r\) 
n 


22n  T(n+1) 


.  2n 

l  erfcrj 


IpM. 

Pn(3) 


(8) 


V  e~^2dg 
0J  (i2nerfc  £)2 


(9) 


For  the  Neumann  problem,  the  integer  value  n  =  0  is  sub¬ 
stituted  into  equation  (9)  giving 


Po(r)> 


Vlr  erf  n 
2  erfcrj 


(10) 


Using  (10)  in  equation  (8) ,  the  Neumann  solution  for  temp¬ 


erature  is 


.  _  erf  Ti 

1  “  erfp 


(11) 
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Evaluation  of  the  integral  (9)  for  a  number  of  integer 
values  2n  -  1,2,..  „  has  been  undertaken  by  the  author  and 
is  presented  in  table  2.1.  To  evaluate  Pn(T)),  the  recur- 

t 

rence  formula 


inerfc  £ 


( in  2erfc  £  -  2^in~^erfc  £ 


is  expanded  term  by  term  and  then  groupings  made  of  the 

-ft  2 

coefficients  of  e  ’  and  erfc  £  in  the  numerator  of  the 


equation 


dPn<?)  =  ,7^ 


_f2 

e  "  d£ 


(i  erfc  £) 


where 


Pn<?> 


u 

V 


and  v 


i2nerfc  £ . 


Hence,  the  numerator  is  considered  to  take  the  form. 


i2nerfc  £ 


i^n'”1erfc  £ 


u 


where  the  necessary  value  of  u  must  be  determined. 


' 

. 
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2n 


Pn(rj) 


0 


pew 


VjT  erf  rj 
2  erfCT] 


1 


2 


3 


4 


5 


6 


_  Vtt~  n 
ierfCTj 


V1^ 

P3/2  ^ 

p2(n) 

P5/2  ^ 

p3(n) 


Vi~^2T‘2  t-i)  _  4^r 

i  erfcrj 

2  ^sTtF"  (2n3  +  3n) 

•  3  * 

l  erfct] 


2  (4n4  +  12ti2  +  3)  _  192 

i  erfcTj 


4  Vtt~  (4ti^  +  20n^  +  15n) 

5 

l  erfcrj 

4  <\Jrr~  (8ri^  +  60n^  +  90t^  +  15) 

l  erfcrj 


60  \[W~ 

i^erfc  0 


TABLE  2.1  TABULATION  OF  P  (rj) 

n 


(f^  n<3 


tv: 


3  ~  r 


Ip)- 


p  r  '9 


* 

I  9 


11 


2o2-2  Quadrature  Solution 

If  2n  is  not  a  positive  integer,  computation  of  the 
temperature  may  be  executed  in  quadrature  form.  Consid¬ 
ering  equation  (4)  and  using  an  integrating  factor,  it 
may  be  written 

j2T)dT|  j2T)dT) 

f  (T))  e  =  4nf  (tj)  e  (12) 

Integrating  (12)  twice  and  imposing  the  boundary  conditions 
(5)  yields  the  solution 


g£fr)  + 

erf|3 


4n 


Q  (t]) 


erf  n 

erfp 


Q(3) 


(13) 


where  Q  ( rj) 


d£d£ . 


o  ^  o  u 

The  Neumann  solution  is  recovered  by  putting  n  =  0. 

Equations  (8)  and  (13)  represent  the  formal  solution 
for  the  temperature  distribution  in  the  ice  region  with 
the  selected  surface  temperature  variation,  9Q(t)  =  Ctn. 

A  computer  programme  in  Fortran  IV  language  was 
written  to  evaluate  the  temperature  function  (13).  The 
computation  was  carried  out  on  an  IBM  7040-1401  machine 
in  the  Department  of  Computing  Science,  University  of 
Alberta.  Appendix  A  provides  the  programme.  The  computer 
data  results  are  plotted  in  figures  2.2,  2.3,  and  2.4. 

In  discussing  the  temperature  function  (13)  ,  it  should 
be  noted  that  Q(tj)  is  bounded  in  the  range 
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FIG.  2.2  THEORETICAL  TEMPERATURE  PROFILES 
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FIG.  2.3 


INTEGRAL  Q(|3)  VERSUS  DEPTH  £ 
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FIG.  2.4  DERIVATIVE  Q  (0)  VERSUS  DEPTH  0 
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0<T)<P 

and  calculations  show  that  the  range  of  the  domain  is 
approximately 

oC  0  <0,765 

where  the  maximum  value  of  (3  =  0,765  corresponds  roughly 
to  a  surface  temperature  of  absolute  zero  or  Gq (t)  =  -492 °F 
in  the  case  of  ice  (refer  to  Appendix  D) .  An  important 
point  to  note  from  Fig.  2.2  is  that  the  temperature  pro¬ 
file  is  nearly  linear  for  small  f3  and  any  value  of  n. 


2 . 3  SOLUTION  FOR  INTERFACE  DEPTH 
2,3-1  Explicit  Solution 

The  dimensionless  ice  depth,  f3,  is  not  known  a  priori 
and  as  a  result,  the  temperature  solution  (8)  is  not  com¬ 
plete  as  presented.  Ice  depth  must  be  determined  from 
an  energy  balance  at  the  ice-water  interface.  In  gen¬ 
eral,  for  the  uni-dimensional  semi-infinite  domain,  the 
balance  is 


(15) 


where  subscripts  1  and  2  represent  the  frozen  and  unfrozen 
regions  respectively.  However,  if  region  2  is  considered 
to  be  at  zero  temperature  initially  and  remains  at  zero 
for  t>o,  the  energy  balance  (15)  reduces  to 
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~K 


dx. 


x=b 


T  db 

PL  dt 


(16) 


where  the  subscripts  are  deleted  for  convenience .  The 
temperature  gradient  at  the  edge  of  the  ice  domain  is 


dx. 


x=b 


eofn  (P) 
2^/ic  t 


(17) 


The  dimensionless  temperature  gradient  at  |3  may  be  ob¬ 
tained  by  differentiating  (8)  which  yields 


fn  (3)  = 


-22nr  (n+l)  i2nerfc  P  p'n(P) 

pnTen 


(18) 


where  Pr  (f3)  is  found  directly  from  equation  (9)  .  Sub¬ 


stitution  of  equation  (18)  into  (17)  and  hence  (16)  gives 
the  ice  front  velocity  as 


db 

dt 


K 


e 


_  _  22n  r(n-f-l)  e  ft 

pL  Pn(P)  i2nerfc  p 


(19) 


where  Pn(p)  may  be  obtained  from  table  2.1. 


If  equation  (19)  is  written  in  the  original  form 


db 

dt 


=  2a  Vic  tn  ■‘•Z2  fn'  (P)  , 


(20) 


where  a  =  -CpC/4L  ,  then  the  dimensional  ice  depth  is 

=  2aVl  fttn-1/2  f'(P)  dt. 

O^ 


b 


(21) 
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Integrating  equation  (21)  by  parts  (and  satisfying  equ¬ 
ation  (4)  at  Tj  “  (3)  gives 


b 


4a  V/c~ 
2n  +  1 


n+1/2 


f  O)  dt 


+  2a 

o 


t2"-^ 


(22) 


f  (P) 


dt 


4a  V/c 
2n  +  1 


^.n+1/2  fn'(p) 


+  2a  I 


O 


Using  to  the  properties  of  ice  on  page  44  ,  and  the 

value  of  C  =  l°F/(hr)n,  calculation  shows  that  a  =-.00085; 

this  immediately  reveals  the  fact  that  the  coefficient  for 

the  third  term  of  equation  (22)  contains  a  smaller  order 
2 

number,  a  .  Therefore,  all  integrals  (including  1^)  con¬ 
taining  the  smaller  coefficient  will  be  neglected.  To 
integrate  1^,  it  is  necessary  to  note  that 


33.  . 

dt 


1 

t3/2 


at"  f '  (p) 


(23) 


which  is  obtained  directly  from  the  definition  of  P  and 
equation  (20) .  Using  this  with 


f '  '  (p)  =  -2  3  f '  (P)  , 


(24) 


an  approximate  series  solution  for  (22)  may  be  found. 
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Dividing  the  result  by  2  V/ct  yields  the  series  solution 
a  =  _2_atnf '  (p)  T  _  p2t~1/2  P2(P2-1)  t"1 

P  2n  +  1  L  n  n (n-1/2)  “  *  *  * 

(25) 

Assuming  the  first  term  is  large  compared  with  all  the 
remaining  terms  (refer  to  Appendix  D) ,  then,  the  solution 
for  dimensions  less  ice  depth  simplifies  to 


P 


2a  tnf n ' (8) 

2n  +  1 


(26) 


Solution  (26)  may  be  obtained  directly  by  integrating  (20) 
with  respect  to  time  and  holding  £  fixed. 

In  the  case  of  the  Neumann  problem,  using  n  =  0 
which  indicates  a  temperature  step  change  at  the  surface, 
equation  (26)  reduces  to 


P 


_o  2 

(  CpC  \  e  ° 

\VFi/  erff3 


(27) 


which  shows  that  the  dimensionless  ice  depth  is  a  constant 
dependent  upon  the  thermal  properties  of  the  frozen  region 
and  the  surface  temperature  step  change.  The  ice  front 
progression  is  therefore  proportional  to  Vt  for  this  special 


case. 
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Observing  the  form  of  equation  (27)  ,  a  useful  dim¬ 
ensionless  form  for  time  is 


n 


4  n 
at 


VtF 


and  therefore  the  explicit  solution  (26)  is  written 


n 


e  =  2(2t\T)  T“  fn 


(28) 


where  f  (p)  is  obtained  from  equation  (18) 


2,3-2  Quadrature  Solution 


If  2n  is  not  a  positive  integer,  the  quadrature 
method  may  be  used  to  calculate  the  ice  depth,  p.  Dif¬ 
ferentiating  (13)  and  evaluating  at  (3  gives 


f  '  o) 

n 

=  4nQ '  (p)  -  ~ 

Vi r 

-p2  r  -| 

erfp  L1  +  4nQ(P)J 

(29) 

I 

R2 

=  e-p  f  fntt) 
oJ 

*2 

where 

Q  (P) 

e"  d£ . 

Substituting  (29)  into  (17)  and  (16) ,  gives  the  ice  front 
velocity 


db 

dt 


-KG 

o 

2pL  *fict 


4nQ ' (p) 


Yir  erfp 


+  4nQ(P) 


(30) 


The  dimensionless  ice  depth,  p,  is  therefore 
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FIG.  2.5  THEORETICAL  PROGRESS  OF  INTERFACE 
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erfp 


( 


1  +  4nQ(g) 

2n  +  1 


T 


n 


(31) 


which  is  found  by  integrating  (30) ,  with  respect  to  time 
holding  p  fixed,  similar  to  the  method  used  in  section  2.3-1. 
The  heat  transfer  rate  at  the  interface  may  be  cal¬ 


culated  by  obtaining  values  of  Q(P)  and  Q  (£)  from  figures 
2.3  and  2.4  respectively  and  substituting  into  equation 
(30)  and  then  (16) . 

For  convenience,  a  plot  of  (3  versus  t  has  been  in¬ 
cluded  as  figure  2.5  in  which  the  convergence  to  the 
limiting  value  of  0.765  is  illustrated.  Although  an 

ice-water  system  has  been  considered  throughout,  figure 
2.5  is  applicable  to  any  medium  changing  phase  with 


9 


o 


2 . 4  SINUS 01 DUAL  TEMPERATURE  AND  PHASE  CHANGE 

A  solution  for  ice  depth  with  a  sinusoidual  surface 
temperature  variation  is  considered  to  be  an  extremely 
important  case.  A  practical  example  is  frost  penetration 
in  soil  or  formation  of  lake  ice  as  governed  by  daily 
(or  annual)  atmospheric  temperature  variations. 


A  power  series  may  be  used  to  mathematically  des¬ 
cribe  the  sinusoidual  surface  temperature  and  it  is  pro¬ 
posed  that  the  temperature  regime  may  be  represented  by 
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oo 


An  tn  fn(rj)  . 


(32) 


n  =0 


The  superposition  of  temperatures  should  provide  a  good 
approximation  to  the  actual  temperature  regime  if  CpC/L 
is  of  small  order.  The  relative  values  of  this  number, 
for  various  materials,  are  shown  in  Appendix  D.  If  the 
latent  heat  of  fusion  is  large  (i.e.  ice) ,  then  the  inter¬ 
face  will  be  relatively  slow  moving,,  If  the  specific 
heat  of  the  solid  is  relatively  small,  then  surface  temp¬ 
erature  variations  will  rapidly  vary  the  temperature 
regime  in  the  solid  region.  Thus  the  superposition  method 
should  be  quite  accurate  for  large  ice  depths  and/or 
for  a  slowly  moving  interface.  In  fact,  using  a  periodic 
boundary  temperature,  experimental  evidence  shows  that 
the  interface  is  virtually  stationary  for  a  major  portion 
of  the  cycle  (Chapter  IV) ;  thus  the  superposition  method 


should  be  applicable  to  this  case.  Considering  the  power 


then,  the  temperature,  in  the  region  0<'n<C£>  is 


(33) 


9(ri) 


where  ©  is  the  maximum  temperature  amplitude.  Equation  (33) 
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satisfies  the  boundary  conditions  (5) .  For  each  half¬ 
cycle,  the  thermal  properties  will  change  as  shown  below: 


Range 

0  <C(JDt  «<  7T 
7T  CDt  2  7T 


Properties 

^1'  cPi'  Pi  ^  ”  freezing  region) 
/c 2  r  CP2'  P2  (2  ”  melting  region) 


For  the  second  half-cycle,  the  boundary  conditions  will 
remain  the  same  and  the  initial  conditions  will  be  nearly 
the  same  as  those  of  the  first  half -cycle.  In  Chapter  IV, 
experimental  evidence  will  be  shown  which  indicates  that 
the  initial  domain  temperature  is  virtually  zero  at  the 
beginning  of  each  cycle  and  for  as  many  as  three  cycles. 
Applying  the  energy  balance  (16) ,  the  ice  front 
velocity  is 


(P)  + 


5 

a 

51 


f5(0)... 


(34) 


where  a  =  ait  and  property  subscripts  1  and  2  are  deleted. 
In  dimensionless  form,  the  ice  depth  is 


-Cp0. 


(3 


m 


2L 


a  f,  (P) 


3 

a 

3  J  7 


(P)  + 


a 


5 ’ll 


fc<0) 


(35) 


Substituting  equation  (18)  appropriately  into  (35)  gives 
the  ice  depth  as 


P 


Cp9  e 
r  m 


-P‘ 


_ a_ _ 

P1(p)i2erfc  p 
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P3  (p) i6erfc  p 
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J  1024 

+  ~TT~ 


5 

a 


P5(P)i10erfc  p 


(36) 


whereby  the  necessary  functions  for  Pn(@)  may  be  obtained 
from  table  2.1. 

For  comparison  of  theory  with  experimental  results, 
only  the  first  three  terms  of  the  series,  vis  equation  (35), 
are  used.  For  consideration  of  the  half-cycle  only,  terms 
beyond  the  third  are  small.  The  first  term  of  the  series 
is  the  solution  for  a  linear  boundary  condition.  In 
applying  equation  (35)  ,  a  more  useful  form  is 


0 


-a  f1  (p) 
3 


3, 
a  f. 


(P) 


3  i  1 


(37) 


where  4>  =  2L/Cp©^.  From  the  computer  data,  a  plot  of  (3 

versus  0  (reciprocal  of  temperature)  has  been  included  as 

figure  2.6.  The  figure  indicates  that  the  maximum  ice 

depth  increases  sharply  for  values  60  or  alternatively, 

* 

as  9^  increases.  An  important  criterion  for  predicting 
maximum  ice  depth  is,  therefore,  0  ,  for  a  mean  temperature 
at  zero  datum  and  for  one  cycle. 


*  A  vast  number  of  field  trials  in  the  study  of  annual  frost 
penetration  have  indicated  that  the  number  of  degree- 
days  below  freezing  per  annum i.e.  proportion  of  area 
of  sinusoidual  curve  below  32  F  as  compared  with  that  above 
32 °F,  will  govern  the  depth  of  freezing.  The  trials  in¬ 
dicated  that  ice  depth  will  increase  if  the  mean  is  30  F 
and  below. 
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FIG.  2.6  RECIPROCAL  TEMPERATURE  <$>  VERSUS  DEPTH  6 


CHAPTER  III 


EXPERIMENTAL  CONSIDERATIONS 


3.1  EXPERIMENTAL  APPARATUS 

During  the  preliminary  design  of  the  apparatus,  the 
primary  consideration  was  to  build  a  device  which  would 
freeze  water  one-dimensionally  from  an  upper  flat  surface 
and  be  able  to  view  the  freezing  interface  laterally 
through  the  side  walls  of  the  device 0  The  original  de¬ 
sign  was  satisfactory  for  this  type  of  ice  formation 
using  a  surface  temperature  starting  initially  at  the 
freezing  point  and  varying  by  a  simple  power  law  of  time0 
Later,  it  was  necessary  to  modify  the  refrigeration  jacket 
by  installing  electrical  heating  coils  to  help  impose  a 
periodic  temperature  variation  at  the  surface  plate.  In 
the  final  phase  of  the  experiments,  a  completely  new  re¬ 
frigeration  jacket  was  designed  and  built  to  undertake 
two-dimensional  ice  formation  tests. 

3.1-1  Original  Test  Cell 

The  diagrammatic  arrangement  of  the  apparatus  is 
shown  in  fig.  3.1.  It  consisted  of  a  transparent-walled, 
rectangular  plastic  tank  (internal  dimensions  10.87  x  7.50 
x  7.94  inches.)  cooled  at  the  top  by  an  81.5  square  inch 
x  1/8  inch  thick  copper  plate  which  formed  an  integral 
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FIG. 3.1  DIAGRAMMATIC  ARRANGEMENT  OF  APPARATUS 
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part  of  a  refrigeration  jacket „  The  refrigeration  jacket 
contained  1/4  inch  O.D.  copper  tubing  bonded  to  the  copper 
plate  by  epoxy  resin*  Freon  12  was  circulated  through 
the  tubing  by  a  vapour-compression  refrigerator.  To  pro¬ 
vide  uniform  heat  transfer  from  the  copper  plate  to  the 
refrigerant,  the  tubing  was  submerged  in  an  ethylene 
glycol/water  solution  capable  of  temperatures  down  to 

O 

-34  F  before  freezing.  The  sides  and  base  of  the  plastic 
(Perspex)  tank  were  insulated  by  3  inch  thick  styrofoam 
with  removable  panels  on  the  front  and  rear  sides  for 

periodic  observation  of  the  ice  growth.  Silicone  rubber 
cushions  were  inserted  at  each  side,  as  shown  in  fig.  3.1, 

to  reduce  the  danger  of  expansion  breakage  in  the  lateral 
direction.  To  prevent . vertical  lifting  of  the  copper  sur¬ 
face,  the  water  was  allowed  to  expand  into  a  stand  pipe. 

The  temperature  of  circulating  refrigerant  in  the 
copper  tubing  was  controlled  by  adjusting  a  pressure¬ 
regulating  valve  in  the  return  line.  The  refrigerant 
temperature  was  held  within  +  1  F  of  the  desired  test 
temperature. 

The  problem  of  differential  contraction  between  the 
copper  plate  and  the  plastic  tank  was  overcome  by  bonding 
a  1/16  x  1/4  inch  hard  rubber  seal  to  the  top  edge  of 
the  tank.  The  copper  surface  plate  rested  on  the  seal. 

A  small  amount  of  highly  water-repellent  silicone  lubricant 
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was  used  on  the  seal.  To  complete  the  assembly,  four 
yellow-brass  threaded  rods  were  installed  near  the  edge 
of  the  plastic  tank  base  and  with  provision  of  a  1/2  inch 
aluminum  plate  cover  on  top  the  refrigeration  jacket,  the 
copper  surface  plate  was  screwed  down  onto  the  seal. 

To  facilitate  removal  of  air  bubbles,  a  1/4  inch  O.D. 
overflow  tube  was  incorporated  in  one  corner  of  the  sur¬ 
face  plate.  Tilting  the  surface  forced  the  air  out  through 
the  tube  during  filling;  the  tube  was  clamped  shut  when 
all  the  air  was  removed.  For  convenience  during  filling 
and  draining,  the  complete  test  cell  was  mounted  over  a 
laboratory  sink  as  shown  in  fig.  3.2. 

3 • 2  TEMPERATURE  INSTRUMENTATION 

Temperatures  were  measured  using  eight  iron-constantan 

0 

thermocouples,  with  a  resolution  of  +  0.25  F  and  calibrated 
at  -40 °F,  0 °F  and  32 °F  with  a  linear  variation  between 
these  points.  They  were  housed  in  Type  316  stainless  steel 
tubes  (0.0625  inch  O.D.  x  0.016  inch  wall)  located  at  discrete 
levels  ranging  between  0.120  and  1.035  inch  below  the 
surface.  The  tubes  were  mounted  parallel  to  the  upper 
surface  (along  isotherms)  to  reduce  conduction  losses; 
they  were  staggered  at  various  levels  in  three  vertical 
rows,  each  row  spaced  2-1/2  inches  from  the  other.  The 
tubes  afforded  good  chemical  and  physical  protection  to  the 
thermocouples  and  prevented  thermocouple  wires  from  sagging. 
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Initially,  there  was  some  question  whether  the 
presence  of  stainless  steel  tubes  inside  the  test  section 
would  distort  the  isotherms  and  create  test  result  errors. 
However,  it  was  considered  that  isotherm  distortion  would 
be  insignificant  when  taking  into  account  the  large  dim¬ 
ensions  of  the  test  cell  as  compared  with  the  small  size 
of  tubing. 

Two  additional  thermocouples  were  installed  in 
small  diameter  holes  drilled  into  the  surface  plate  from 
the  upper  side.  Temperature  measurements  from  these 
were  averaged  to  furnish  the  experimental  "surface" 
temperature. 

During  preliminary  tests,  two  strip  chart  recorders 
were  used  to  record  thermocouple  temperatures.  A  high 
sensitivity  single  channel  Honeywell  Electronik  19  re¬ 
corder  was  connected  to  one  surface  thermocouple  to 
evaluate  the  maximum  "oscillation"  amplitude  of  surface 
temperature.  This  was  necessary  because  the  time  de¬ 
pendent  temperature  boundary  conditions  were  obtained  by 
approximating  the  desired  function  of  time  with  a  series  of 
small  step  changes.  For  larger  step  changes,  the  oscil¬ 
lation  was  greater.  The  maximum  amplitude  of  this  oscil¬ 
lation  was  of  the  same  order  as  the  refrigerant  tempera¬ 
ture  variation.  A  lower  sensitivity  sixteen  channel 

o  o 

Brown -Honeywell  recorder,  with  a  range  of  -50  F  to  +100  F 
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provided  a  chart  record  of  the  remaining  nine  thermo¬ 
couples  during  each  run,  the  circuit  of  sixteen  channels 
being  recorded  every  2-1/2  minutes . 

3  o  3  INTERFACE  DEPTH  MEASUREMENT 

To  facilitate  ice  depth  measurement,  a  flat  aluminum 
bar  was  mounted  laterally  across  the  front  of  the  clear 
plastic  test  section,,  Threaded  brass  rods  were  installed 
at  each  end  of  the  bar  which  allowed  uniform  vertical 
movement  by  rotating  each  rod  an  equal  number  of  turns . 

Each  turn  moved  the  bar  vertically  CL  020  inch.  Scribed 
micrometer  knobs  and  revolution  counters  were  fitted  to 
the  rods.  A  rifle-sight,  with  fore  and  aft  vertical  ad¬ 
justments,  was  mounted  in  a  lateral  slot  on  top  of  the 
aluminum  bar.  The  rifle-sight  was  zeroed  on  the  under¬ 
side  of  the  surface  plate  prior  to  each  test  and  after 
repeated  checks  the  accuracy  of  sighting  on  the  surface 
was  found  to  be  +  0.002  inch.  However,  the  overall 
accuracy  of  the  arrangement  was  considered  to  be  +  0.005  inch. 


3.4  MODIFICATIONS  TO  EXPERIMENTAL  APPARATUS 
3.4-1  Slab  Test 

A  simple  addition  was  made  to  the  test  cell  to  simu¬ 
late  freezing  of  a  shallow  lake.  The  requirement  was  to 
provide  a  finite  water  layer  over  a  solid  substrate  and 
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to  use  a  power  law  variation  of  surface  temperature. 
The  choice  of  material  for  the  substrate  was  styrofoam 
which  provided  a  near  zero  heat  transfer  boundary  con¬ 
dition,  The  styrofoam  properties  were  (approximately) ; 


0,02  BTU/hr  f t °P 

0,2  BTU/lb  °F 

1,9  

2 

0,0525  ft  /hr  (calculated) 


K 


C 


P 

P 


/c 


The  substrate  was  pressed  firmly  into  the  lower  section 
of  the  plastic  tank  and  the  remaining  volume  was  filled 
with  tap  water  to  a  depth  of  1  inch  in  preparation  for 
the  test. 


3.4-2  Sinusoidual  Temperature  Tests 

After  the  initial  set  of  tests,  using  the  simple  power 
law  variation  in  surface  temperature,  a  periodic  tempera¬ 
ture  variation  over  three  cycles  was  attempted.  External 
heat  lamps  were  required  to  raise  the  surface  tempera¬ 
ture  to  the  necessary  degree  above  freezing  prior  to  re¬ 
turning  to  below  freezing  in  each  cycle.  For  a  second 
test,  an  electrical  heating  wire,  wound  over  two  3/4  inch 
O „ D.  Teflon  rods,  was  submerged  in  the  glycol  solution 
in  the  refrigeration  jacket.  This  heater  replaced  the 
cumbersome  heat  lamp  arrangement.  Electrical  power  to  the 
heating  element  was  controlled  by  a  Variac  unit. 
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3.4-3  Two-Dimensional  Tests 

With  completion  of  the  one-dimensional  tests,  a 
revised  refrigeration  jacket  was  manufactured  and  mounted 
onto  the  original  plastic  tank  for  a  study  of  two  dim¬ 
ensional  ice  formation..  For  this  particular  study,  ex¬ 
tensive  modifications  were  required.  The  complete  two- 
dimensional  test  apparatus  is  shown  in  fig.  3.2. 

The  aim  of  the  test  was  to  locate  the  steady  state 
ice  profile  under  the  edge  of  a  surface  held  at  a  uni¬ 
form  temperature  below  freezing  and  with  the  adjacent 
surface  held  at  a  uniform  temperature  above  freezing. 

The  revised  jacket  consisted  of  two  copper  plate 
surfaces  separated  by  an  insulation  of  1/4  inch  asbestos 
board.  The  yellow-brass  side  walls  of  the  jacket  were 
identical  to  those  used  In  the  original  test  cell.  Id¬ 
entical  copper  cooling  coils  were  mounted  on  each  sur¬ 
face  plate  and  connected  in  parallel  to  the  refrigerator 
with  a  shut-off  valve  incorporated  to  control  Freon  12 
flow  to  the  right-half  coil.  An  important  additional 
feature  was  the  inlet/outlet  ports  in  the  right  half  of 
the  jacket  for  circulation  of  an  externally  heated  glycol/ 
water  solution  for  maintaining  a  uniform  temperature 
above  freezing. 

Five  iron-constantan  thermocouples  were  installed 
in  the  jacket,  two  on  the  cooled  surface,  one  centered 
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FIG.  3.2  VIEW  OF  APPARATUS  FOR 

TWO-DIMENSIONAL  EXPERIMENTS 
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in  the  asbestos  board,  and  two  on  the  heated  surface,, 

3  o  5  EXPERIMENTAL  METHOD 

The  analytic  model  for  the  power  law  variation  outlined 
in  Chapter  II  was  applicable  only  to  the  case  when  the 
water  was  initially  at  a  temperature  of  32 °F.  Consequently, 
prior  to  each  test  run,  the  water  was  cooled  from  room 
temperature  to  a  point  when  the  temperature  profile  was 
steady  and  as  close  to  the  zero  datum  as  possible.  This 
required  approximately  3  hours ,  at  which  time  the  lower 
thermocouple  in  the  water  registered  a  steady  35-36  °F 
which  was  considered  to  be  the  limiting  case  for  the 
apparatus  with  only  3  inches  of  styrofoam  insulation.  The 
initial  departure  from  zero  datum  as  described  above  was 
considered  as  a  condition  of  "superheat". 

3o5-l  Temporal  Origin  and  Supercooling 

During  the  initial  power  law  tests ,  it  was  discovered 
that  pronounced  supercooling  of  the  water  occurred.  This 
was  a  departure  from  the  theoretical  requirement  for  an 
initial  temperature  at  zero  datum.  The  difficulty  was 
surmounted,  in  subsequent  power  law  tests,  by  producing 
a  thin  layer  of  ice  on  the  surface  and  then  allowing  thaw- 
back  until  only  an  extremely  thin  film  of  ice  was  visible 
on  the  surface.  The  test  run  proper  was  then  started  from 
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Even  with  precautions  to  overcome  the  supercooling 
effect,  the  temporal  origin  of  each  test  run  was  difficult 
to  fix*  To  provide  a  uniform  basis  for  comparison  of 
experiment  with  the  analytic  model,  the  temporal  origin 
was  considered  to  be  a  time  point  found  by  rearward  ex¬ 
trapolation  of  the  surface  test  temperatures  to  32 °F. 

Since  thermocouple  temperatures  were  measured  over  a 
2-1/2  minute  cycle  and  all  tests  were  more  than  one  hour 
duration,  the  extrapolation  was  sufficiently  accurate 
to  establish  the  temporal  origin* 

3 o5-2  Typical  Test  Procedure 

Prior  to  each  test,  additional  water  was  forced  into 
the  enclosed  test  cell  to  remove  air  bubbles  which  accumu¬ 
lated  at  the  surface  plate  between  tests*  The  bulk  of 
the  water  was  then  cooled  for  several  hours  to  near  32  °F* 
Following  this,  a  thin  ice  layer  was  produced  by  quickly 
reducing  the  surface  temperature  below  the  freezing  point* 

O 

The  surface  temperature  was  then  reset  to  32  F  and  thaw- 
back  allowed  as  described  previously.  During  thawback , 
the  ice  depth  measuring  device  was  "zeroed"  on  the  sur¬ 
face  plate  by  removing  the  insulated  observation  panels. 
The  temperature  recorders  were  turned  on  and  at  the 
"critical"  moment  of  thawback,  the  desired  surface  temp¬ 
erature  variation  was  initiated*  Temperature  measurements 
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were  continuous  and  automatic  throughout  test  runs  whereas 
ice  depth  measurements  were  on  a  periodic  inspection  basis. 

3.6  SUPPLEMENTARY  SUPERCOOLING  EXPERIMENTS 

Another  set  of  experiments  were  undertaken  to  find  a 
relation  between  the  surface  cooling  rate  and  supercooled 
time  period  required  for  ice  to  form.  As  in  previous 
experiments ,  tap  water  was  used.  The  apparatus  was  a 
3M  Brand  Thermoelectric  Cooler  Model  11  B?  it  consisted 
of  a  0.32  cubic  inch  cylindrical  cavity  insulated  with 
styrofoam  and  cooled  by  thermoelectric  modules.  Thermo¬ 
couples  and  an  electrical  power  source  were  installed 
for  the  tests.  By  comparison,  the  volume  of  the  freezing 
cell  used  in  the  main  experiments  was  over  2,000  times 
larger  than  the  thermoelectric  cell. 

A  set  of  10  test  runs  were  conducted  with  the 
thermoelectric  unite  Results  were  compared  with  three 
prior  tests  on  the  large  freezing  cell. 
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CHAPTER  IV 


DISCUSSION  OF  RESULTS 


4. 1  SUPERCOOLING 

One  of  the  first  problems  encountered  in  the  experi¬ 
mental  programme  was  that  of  water  supercooling  near  the 
surface  prior  to  formation  of  ice„  During  initial  tests, 
water  temperatures  as  low  as  22  °F  were  recorded  before 
crystallization..  Numerous  ice  spikes  or  "spicules"  pro- 

v 

pagated  rapidly  downward  into  the  water?  in  one  test, 
spicules  penetrated  3/8  inch  within  a  period  of  30  seconds. 
Following  their  formation,  the  lower  ends  of  the  spicules 
melted  back  uniformly  and  gradually  degenerated  into  a 
thinner  but  homogeneous  layer  of  ice.  From  this  point, 
a  plane  ice-water  interface  propagated. 

A  typical  temperature  record  (n  -  1.0  and  C  =  25.0°F/hr.) 
indicating  supercooling  is  shown  in  fig.  4.1  and  clearly 
shows  a  metastable  period  followed  by  a  very  rapid  sur¬ 
face  temperature  increase.  This  is  compared  with  another 
test  (n  “  0.71  and  C  -  25.0°F/hr*  ±)  where  ice  was  formed 
initially  and  melted  back  to  a  thin  film  prior  to  the 
test  run.  Fig.  4„2  illustrates  the  steps  (a)  to  (d)  which 
are  believed  to  occur  during  supercooling  and  subsequent 
ice  formation. 

A  second  test  (n  =  1.0  and  C  =  25.0°F/hr)  was  carried 
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FIG. 4. 1  TYPICAL  SURFACE  TEMPERATURE  RECORDS 
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FIG. 4. 2  TEMPERATURE  PROFILES  DURING  SUPERCOOLING 
AND  ICE  FORMATION 
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out  but  with  supercooling  eliminated  as  previously  des¬ 
cribed,,  Within  the  limits  of  experimental  error  and  one 
hour  after  the  time  zero  as  defined  in  section  3.5-1,  no 
appreciable  difference  could  be  measured  in  ice  depth  with 
and  without  supercooling,,  This  leads  to  the  important 
point  that  the  temporal  origin  in  experiments  should  not 
be  taken  as  the  time  when  ice  just  forms  but  rather  as 
the  time  defined  in  section  3.5-1.  For  longer  time  per¬ 
iods,  after  supercooling  with  subsequent  ice  formation, 
the  analytic  solution  should  provide  an  accurate  predic¬ 
tion  of  ice  depth.  Fig.  4.1  further  substantiates  the 
idea  that  supercooling  does  not  cause  a  serious  departure 
from  the  analytic  model?  for  longer  times,  the  actual 
surface  temperature  with  supercooling  converged  to  the 
nominal  temperature. 


4.1-1  Results  of  Supplementary  Supercooliy  Experiments 
For  the  particular  volume  of  test  cavity  in  the 
thermoelectric  unit,  a  definite  relation  was  found  be¬ 
tween  the  surface  temperature  cooling  rate,  d©o/dt,  and 
the  maximum  time,  t  ,  of  supercooling  of  water  prior  to 
crystallization.  The  empirical  relation  was  found  to  be 
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Considering  supercooling  and  the  particular  linear  sur¬ 
face  temperature  shown  in  fig.  4.1,  the  empirical  relation 
predicts  that  tc  -  0.304  hr.  using  0q  ^  Ctn  where  n  =  1.0; 
this  agrees  well  with  the  time  of  crystallization  in  fig. 
4.1.  In  more  general  terms,  the  maximum  time  (using  the 
above  surface  temperature)  is 


nC 

1.22 


1 

3 . 54-n 


which  points  out  the  possibility  of  no  ice  forming  when 

n  -  3.54  and  thus  t  oo . 

c 

During  the  small  cavity  tests,  the  sudden  temperature 
rise  associated  with  formation  of  stable  ice  was  constant 
at  12.75  +  0.5  °F  for  all  test  runs  regardless  of  the  de¬ 
gree  of  supercooling.  This  was  contrasted  with  a  typical 

o  t 

surface  temperature  rise  of  5  to  7  F  in  the  larger  appar¬ 
atus.  The  difference  stems  from  the  difference  in  dim¬ 
ensions  of  the  two  systems  and  the  volume  of  ice  formed. 
Initially,  it  was  thought  that  stable  ice  filled  the 
smaller  unit  almost  instantly  with  no  apparent  thawback. 
However,  the  latent  heat  released  was  calculated  to  be 

only  .015  BTU  whereas  to  fill  the  unit  completely  with 

. 

ice  would  require  a  release  of  0.171  BTU.  A  period  of 
up  to  3  minutes  was  required  for  stable  ice  to  form,  after 
thawback,  in  the  larger  unit. 

The  degree  of  supercooling  encountered  in  the 
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supplementary  experiments  ranged  as  low  as  21 °F  below 
the  normal  freezing  temperature  using  ordinary  tap 


water o  Turnbull 


13 


has  indicated  that  under  special 


conditions ,  very  pure  water  can  be  supercooled  to  as 
low  as  -40 °F  without  crystallizing  into  ice. 


40 2  INITIAL  ICE  FORMATION 

From  the  beginning  of  the  test  programme,  when 
difficulty  was  encountered  in  trying  to  establish  a 
criterion  for  temporal  origin,  a  routine  was  formed 
to  observe  the  initial  ice  formation.  These  observa¬ 
tions  revealed  that  a  small  speck  of  ice  (possibly 
several  actual  crystals)  would  form  on  the  surface 
plate  at  a  random  point.  Shortly  thereafter,  similar 
specks  would  form  at  other  points.  Almost  immediately 
a  crystal  speck  formed,  additional  crystals  would 
attach  themselves  to  it  fanning  out  quickly  only  in  a 
radial  horizontal  direction  along  the  plate  surface. 
There  appeared  to  be  no  appreciable  crystal  growth  in 
the  vertical  direction  until  the  complete  plate  sur¬ 
face  was  covered  with  an  initial  thin  film  of  ice 
crystals  * 
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40  3  RESULTS  USING  POWER  LAW  SURFACE  TEMPERATURE 
4  o  3 -1  Ice  Depth 

Visual  studies  of  the  ice-water  interface  showed  it 
to  be  planar,  thus  implying  that  the  ice  crystal  growth 
on  the  interface  proceeded  uniformly  in  a  vertical  dir¬ 
ection  o  Within  the  limits  of  accuracy  of  the  measuring 
device,  it  was  estimated  that  departure  from  flatness 
was  of  the  order  of  0*01  inches*  It  is  uncertain  as 
to  how  the  growth  proceeded  immediately  adjacent  to 
the  transparent  walls*  Removing  the  observation  panels 
during  depth  measurement  periods  caused  thawback  in  this 
region  but  allowed  accurate  determination  of  the  ice 
depth  over  the  main  interface*  The  estimated  thawback 


adjacent  to  the  walls  was  usually  about  0*1  inch  but 
this  was  small  compared  with  the  total  interface  area 
of  81.5  square  inches.  This  suggests  that  one-dimen¬ 
sional  studies  could  be  executed  readily  with  a  much 
smaller  surface  area  and  volume  of  water. 

Theoretical  and  experimental  results  for  the  varia¬ 
tion  of  interface  depth  with  time  are  shown  in  fig.  4.3 
in  dimensional  form.  They  are  calculated  on  the  assump¬ 
tion  that  the  thermal  properties  of  ice  (over  the  temp¬ 
erature  range  of  -20  °F  to  +32' °F)  are  constant  at: 


K  =  1.4  BTU/hr  f  t  °F 
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FIG.  4. 3  INTERFACE  DEPTH  VERSUS  TIME  (POWER  LAW) 
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Fig0  4.3  illustrates  the  effect  of  varying  the  bound¬ 
ary  temperature  index  and  indicates  generally  good  agree¬ 
ment  between  theory  and  experiment.  It  is  interesting 
to  note  the  large  departure  from  the  Neumann  Solution 
(n  “  0)  indicating  its  inadequacy  for  predicting  ice 
depth  in  circumstances  where  n  /  0. 

4,3-2  Temperature  Profile  In  Ice 

Typical  theoretical  and  experimental  temperature 
profiles,  within  the  ice  layer,  are  shown  in  fig,  4,4 
for  n  ~  0,71  and  C  =  25. 0  °F/hr^°  The  value  of  the 

constant,  C,  may  be  found  quickly  for  any  power  law 
variation  since  it  is  numerically  equal  to  the  tempera- 

O 

ture  reduction  below  32  F  (zero  datum)  when  the  time, 
t,  is  equal  to  1  hour.  Good  agreement  with  the  theoret¬ 
ical  temperature  profile  is  evident.  Any  discrepancy 
between  the  theoretical  and  experimental  ice  depth  may 
be  attributed  to  three  factors s  (1)  the  choice  of  temp¬ 
oral  origin,  (2)  the  uncertainty  in  the  choice  of  thermal 
properties  (especially  considering  a  40 °F  temperature 
variation  within  the  ice  layer)  and,  (3)  the  degree  of 
superheat  within  the  liquid  domain.  Fig.  4.4  confirms 
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FIG .  4o  5  VARIATION  OF  TEMPERATURE  WITH  DEPTH 
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that  the  temperature  profile  is  almost  linear  within 
the  ice. 

4.3-3  Slab  Experiment 

A  slab  experiment  was  undertaken  with  1"  depth  of 
water  over  a  styrofoam  substrate.  Comparison  is  made 
with  theory,  for  n  =  1.0  and  C  =  27.0°F/hr,  only  up  to 
thetime  when  the  water  is  completely  frozen  to  the  top 
of  the  substrate.  After  that  time  the  problem  becomes 
one  of  pure  conduction  in  a  stratified  medium  with  an 
initial  non-zero  temperature  distribution  in  the  upper 
stratum.  Fig.  4.5  clearly  reveals  the  effect  of  an  in¬ 
sulating  substrate  such  as  could  exist  beneath  a  shallow 

* 

lake;  it  shows  that  the  interface  accelerates  upon  approach¬ 
ing  the  insulated  substrate  and  also,  that  the  theoretical 
prediction  is  satisfactory  until  3/4  of  the  water  is  frozen. 
At  this  time,  the  absence  of  a  latent  heat  "sink"  causes 
a  pronounced  change  in  the  rate  at  which  the  bulk  ice 
temperature  decreases.  Unfortunately,  the  substrate  was 
not  instrumented  for  temperature  and  hence,  temperatures 
are  not  available  for  that  region. 


4.4  RESULTS  USING  PERIODIC  SURFACE  TEMPERATURE 

Two  experiments  were  tackled  using  a  periodic  surface 
temperature.  The  approximate  test  conditions,  averaged 


*  Refer  to  Lightfoot 
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over  three  complete  cycles,  weres 

Temperature  Temperature  od 

Test  NOo  Amplitudes  Mean  Radians/hr 

1  +  8°F  +  34  °F  t r 

2  +  19  F  +  33  F  t r 

For  theoretical  comparison,  fig.  2.6  is  used?  it  is  based 
on  a  sinusoidual  temperature  symmetrical  about  a  mean 

O 

of  +32  Fa  Fig*  4.6  compares  the  actual  surface  tempera¬ 
ture  with  the  desired  sinusoidual  temperature  for  test 
no.  2.  Test  no.  2  was  considered  more  useful  for  com¬ 
parison  because  the  ice  depth  measurements  were  more 
accurately  done  than  those  of  test  no.  1. 

4.4-1  Ice  Depth  For  Periodic  Case 

Following  the  first  complete  cycle  of  freezing,  a 
thin  layer  of  ice  remained  below  the  surface.  However, 
during  test  no.  1,  no  provision  was  made  in  the  tank  to 
hold  the  layer  at  the  level  where  it  formed  and  consequent¬ 
ly,  it  disintegrated  easily  after  melting  back  to  1/16 
inch  thickness?  this  resulted  in  broken  ice  sheets  float¬ 
ing  to  the  top  of  the  tank.  For  test  no.  2,  the  problem 
was  overcome  by  installing  four  1/8  inch  threaded,  stain¬ 
less  steel  rods  mounted  vertically  inside  the  test  tank 
to  provide  anchorage  for  the  ice  layer  formed.  The  shape 
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FIG.  4.6  SURFACE  TEMPERATURE  FOR  PERIODIC  CASE 
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of  the  ice  formation  near  each  rod  indicated  that  the 
thermal  regime  was  influenced  1/8  inch  radially  outward 
from  the  rod  surface »  Hence  the  presence  of  the  rods 
did  not  greatly  affect  the  overall  ice  depth„ 

For  a  selected  sinusoidual  temperature  amplitude 

O 

of  +  18  F,  the  reciprocal  freezing  temperature  is 
$  =  32 o 8.  Fig.  4.7  shows  the  experimental  ice  depth 
(test  no.  2) ,  clearly  illustrating  the  regions  of  freez¬ 
ing  and  melting,  over  three  complete  cycles.  Using 
these  results  and  fig.  2.6,  a  comparison  of  "theory" 
versus  experiment  is  given  in  table  4.1. 


SINUSOIDUAL  CYCLE  DATA _ VALUE  OF  3  (ICE  DEPTH) 


cot“a 

(Radians) 

Sin  a 

6o 
(  F) 

Theory  Cycle 

1  Cycle  2 

Cycle  3 

0.2 

.199 

3.6 

.045 

.000 

.055 

.065 

0.4 

.389 

7.0 

.063 

.027 

.082 

.094 

0.6 

.565 

10.2 

.077 

.065 

.096 

.109 

0.8 

.717 

12.9 

.088 

.090 

.105 

.118 

1.0 

.842 

15.2 

.096 

.106 

.111 

.122 

1.2 

.932 

16.8 

.104 

.114 

.115 

^24 

1.4 

.985 

17.7 

.110 

.119 

a) 

0) 

n  h 

1.6 

.999 

18.0 

.116 

.123 

T3  rH 
<0 

(0 

(0 

1.8 

.974 

17.5 

.119 

.126 

a) 

U  H 

cn  a> 

JH  rH 

Q)  V 

<u  u 

2.0 

.909 

16.4 

.123 

.127 

>i  >i 

m  u 

>i  >i 

to  u  » 

2.4 

.674 

12.1 

.128 

.127 

rH 

a 

rH  CM 

a 

Q)  o, 

o)  o  rf 

3.0 

.140 

2.5 

.132 

.120 

u  u 
h  m 

U  >H  G 

h  m  fd 

r 

TABLE  4.1  "THEORETICAL"  VERSUS  EXPERIMENTAL  DEPTH 

i 

OF  FREEZING  FOR  PERIODIC  CASE 


£5 


i&rtt  bs>3BO.lbal  fcoi  no**  rroiit>DQ'Soi  sol  ©rto  So 

s ..  ■'*'••  T )  \  :\  9  Ip  i ..  near  r 


.rfyqsl  ©•>!  lit  tsvo  »<£j  *?*£>•  3on  bib 


1  i  :•  '  Sv  i  >c-  ::•  If.  t  :  </  ’.  :  ;> 

J*  -li  -K(X  5  ■  >f  ?  '  >  fc»jL  l  ■'■• 


V  si* .!  -•  ■  "  ..  e  w  , 


£20. 

<-  30. 


CXI. 

CSX  ‘ti 


eei. 


eee. 

*.0 

S£3  , 

see. 

eoe. 


xmaa  jArtfSMweraxa  su^xv  "  jADi^aaosHT  *  i  a  sjsat 


u  “r»  Ti  .  ;  -3  q  j  crc  >  '  > 


53 


The  location  of  the  melting  interface  for  the  second 
half-cycle  is  calculated  by  using  the  properties  for  water 
in  the  evaluation  of  0.  Typical  water  properties  are: 

C  «  1. 005  BTU/lb  °F 

P  m 

L  =  144  BTU/lbm 

/c  «  0o  0051  ft2/hr 

Hence,  for  an  amplitude,  9^  =  +  18 °F,  the  reciprocal 
melting  temperature  is  0  =  15c 9,  Calculation  of  dimen¬ 
sionless  ice  depth  is  by 


p  ice 


b 

2V  Kjt 


2.23b 

Aft 


whereas  calculation  of  dimensionless  melt  depth  is  by 


p  water 


b 

2YV 


7. 00b 


Table  402  furnishes  the  results  for  depth  of  melting. 
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(Based  on  Time  Point) 


SINUSOIDUAL 

CYCLE 

EXP. 

MELT  DEPTH 

THEORY 

EXPERIMENT 

a=oit 

t 

9 

o 

bi 

b2 

b3 

13 

p3 

rad 

hr 

F 

ft 

ft 

ft 

0.2 

.064 

3.6 

— 

— 

— 

.  066 

— 

— 

— 

o 

0 

4^ 

.127 

7.0 

.  000 

.002 

.002 

.091 

.000 

.039 

.039 

O 

0 

CT> 

.191 

10.2 

.000 

.008 

.007 

.110 

.000 

.128 

.112 

05 

0 

o 

.255 

12.9 

.004 

.011 

.011 

.125 

.  055 

.152 

.152 

H* 

0 

o 

.318 

15.2 

.009 

.014 

.015 

.138 

.112 

.174 

.186 

1.4 

.445 

17.7 

.015 

.019 

.020 

.158 

.157 

.199 

.209 

1.6 

.510 

18.0 

.018 

.022 

.023 

.165 

.176 

.215 

.225 

to 

• 

o 

.636 

16.4 

.023 

.025 

.027 

.176 

.202 

.219 

.236 

2.4 

.765 

12.1 

.026 

.028 

.030 

.181 

.208 

.224 

.239 

3.0 

.955 

2.5 

.031 

.032 

.033 

.184 

.222 

.229 

.236 

TABLE  40  2  "THEORETICAL"  VERSUS  EXPERIMENTAL  DEPTH  OF  MELTING 

FOR  PERIODIC  CASE 

Referring  to  fig.  4.7  and  tables  4.1  and  4.2,  the 
depth  of  freezing  will  be  greater  than  the  depth  of  melt¬ 
ing?  this  arises  from  the  difference  in  the  averaged  thermal 
properties  of  ice  and  water.  Generally  this  difference 
exists  for  most  materials.  The  difference  in  freezing  and 

melting  depths  will  occur  using  a  sinusoidual  surface 

* 

temperature  symmetrical  about  a  mean  temperature  at  zero 
datum.  Theoretically,  equation  (34)  indicates  that  for 
assumed  constant  density  across  the  interface,  the  ratio 
of  freezing  depth  to  melting  depth  is  proportional  to 
K yty/c  ignoring  changes  in  the  initial  conditions.  Using 
the  thermal  properties  listed  previously,  the  proportion 


*  Refer  to  page  24  footnote. 
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is  calculated  as  1,38s 1,00,  Considering  the  first  cycle 
shown  in  £ig„  4,7,  the  measured  ratio  is  1,59s 1,00?  the 
discrepancy  may  be  due  to  imperfect  test  conditions  and 
variable  density  across  the  interface. 

The  departure  from  the  analytic  model  will  be  pro¬ 
gressively  greater  as  ice  thickness  increases  because 
the  ‘“initial"  conditions  will  change  slightly  for  each 
cycle.  For  large  depths,  the  ice  bulk  will  behave  as 
a  "semi -infinite"  domain  with  an  "active"  freezing/melt¬ 
ing  zone  near  the  surface, 

4,4-2  Temperature  Profiles  For  the  Periodic  Case 

It  was  observed  that  during  the  second  and  third 
cycles,  no  time  delay  (supercooling)  existed  prior  to 
formation  of  a  new  ice  layer.  In  the  first  cycle,  no 
ice  was  present  initially,  whereas,  in  subsequent  cycles, 
an  initial  ice  layer  was  present  within  approximately  1 
inch  of  the  surface.  This  tends  to  substantiate  the  theory 
that  no  supercooling  will  occur  if  a  few  stable  ice 
crystals  are  available  somewhere  close  in  the  domain. 

An  analytic  consideration  was  that  the  initial  temp¬ 
erature  in  the  domain  (ice  or  water)  must  be  at  32 °F. 

Table  4,3  clearly  shows  that  this  was  a  valid  assumption. 
Experimental  temperature  profiles  are  shown  for  the  third 
complete  cycle  of  tests  no.  1  and  2,  Initial  temperature 
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conditions  are  detailed  at  the  zero  and  one  hour  time 
points o  At  zero  time  in  each  test,  the  position  of  the 
suspended  ice  layer  remaining  from  the  first  two  cycles 
may  be  noted  in  the  temperature  profile0  The  elevations 
of  these  layers  were  (at  time  zero) s 

Test  no,  1  -  ice  layer  .021  to  .024  ft.  below 

surface 

Test  no,  2  -  ice  layer  .032  to  .060  ft.  below 

surface 

The  temperature  oscillation  at  zero  time  of  test  no.  1 
is  not  understood.  With  the  experience  from  test  no.  1, 
the  second  experiment  was  more  carefully  undertaken  and 
hence,  experimental  justification  of  the  initial  condition 
is  based  primarily  on  test  no.  2. 
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CJDt  , 

APPROX.  SURFACE 

THERMOCOUPLE  NO. 

DIME 

TEMP.  3 

4 

5 

6 

7 

8 

9 

10 

Remark 

0 

0 

32.0 

32.2 

36.2 

32.2 

33.7 

44.2 

32.5 

35.0 

45.0 

Test 

w/2 

h 

26.0 

30.2 

33.5 

32.2 

34.5 

41.0 

32.5 

36.2 

43.2 

No.  1 

W 

l 

32  o  0 

32.0 

32.0 

32.0 

32.5 

37.0 

32.0 

33.5 

39.0 

3tt/2 

i-h 

41.7 

32.2 

32.2 

32.0 

32.5 

37.2 

33.0 

/it 

33.0 

39.5 

0 

0 

32.0 

32.5 

33.0 

33.5 

32.0 

32.0 

> 

•H 

33.0 

34.0 

Test 

w/2 

H 

14.0 

17.0 

20.0 

21.0 

25.5 

29.5 

•P 

<0 

i_i 

32.5 

33.0 

No.  2 

w 

i  . 

32.0 

31.5 

32.0 

31.5 

31.5 

31.5 

<u 

a 

32.0 

32.5 

3ir/2 

i-h 

52.0 

39.0 

33.0 

33.0 

32.0 

32.0 

0 

G 

H 

32.0 

33.0 

Thermocouple 

Depth 

.0100 

.0196 

.0304 

.0404 

.0525 

.0629 

.0738 

.0863 

ft. 

TABLE  40  3  TEMPERATURE  PROFILES  FOR  PERIODIC  CASE 
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In  the  two-dimensional  experiments,  two  adjacent 
surfaces  were  maintained  at  a  uniform  temperature?  one 
surface  was  held  at  a  temperature  above  freezing  and  the 
other  below  freezing.  Only  three  experiments  were  under 
taken.  Firstly,  for  a  particular  temperature  difference 
the  ice  profile  was  allowed  to  form  from  an  initially 
melted  domain  and  then, it  was  compared  with  the  profile 
obtained  after  thawback  from  an  initially  frozen  domain 
primarily  to  establish  that  the  profile  had  reached  a 
stable  curve.  Secondly,  the  surface  temperatures  were 
lowered  to  develop  another  ice  profile  to  illustrate 
any  trend  over  the  temperature  range.  Combined  with 
this  category  was  the  comparison  with  a  graphical  solu¬ 


tion,  Brown 


14 


,  which  is  used  to  establish  the  32  F 


isothermal  with  no  regard  to  phase  change  or  variations 
in  properties.  Finally,  from  temperature  measurements 
at  discrete  points ,  the  approximate  location  of  iso¬ 
therms  were  obtained. 


4.5-1  Freezing  Versus  Thawback  Ice  Profile 

The  co-ordinates  of  the  interface  were  measured 
easily  by  sighting  along  the  plane  of  the  interface. 
At  every  co-ordinate,  the  interface  was  exceptionally 
straight  along  the  third  axis.  For  reference,  the 
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temperature  of  the  surface  below  freezing  is  denoted  by 
and  the  surface  temperature  above  freezing  by  T^o 
The  temperatures  in  the  thawback  and  freezing  tests  were? 

(a)  thawback  from  a  frozen  domain;  =  1905°F, 

T2  =  59.  5  °F 

(b)  freezing,  from  a  melted  domain  =  20.5°F, 

T2  =  55 . 5  °F. 

Fig0  40 8  shows  the  resulting  ice  curves?  nearly  identical 
profiles  were  obtained  for  the  given  temperature  differ¬ 
ence,  whether  thawing  back  or  freezingQ  Hence,  the  ice 
profile  was  stable  for  both  tests.  The  plot  also  shows 
that  the  ice  did  not  reach  the  centre  of  the  1/4  inch 
wide  asbestos  insulation  board.  This  is  a  result  of 
the  linear  temperature  increase  across  the  thickness 
of  the  insulation?  a  sharp  temperature  step  change  would 
be  achieved  only  with  a  perfect  thermal  barrier  of  zero 
thickness . 

4.5-2  Two-Dimensional  Ice  Profiles 

Fig.  4.9  shows  the  variation  in  the  ice  depth  with 
a  change  in  surface  temperatures.  The  thinner  ice  layer 
is  a  replot  from  fig.  4.8,  where  T^  =  19.5  F  and  T2  —  59.5  F. 
The  thicker  ice  layer  resulted  from  a  surface  temperature 
of  T1  -  10 . 5  °F  and  T2  =  36.5°F.  Superimposed  on  the 
profiles  in  fig.  4.9  are  typical  32 °F  isotherms  found 


. 


.•  xeqr;.'- j  p  '  >  ;o'v  nit  >  v. 


.  ?  ff-j-  ."  '>  '  ad  ■  '  fi 


' 


t&x  j.i  v  F  *vz-:.y.  3  ?••  •  • 


62 


by  the  graphical  method,.  These  isotherms  would  result 
in  a  single  phase  (steady)  system?  they  are  based  on  a 
"geothermal"  temperature  gradient  of  4 °F  per  inch  which 
is  close  to  the  temperature  gradient  observed  initially 
in  the  water.  Even  though  the  basis  for  the  isothermal 
plot  may  not  be  truly  representative  of  the  particular 
system  (for  a  single  phase) ,  it  does  reveal  that  the 
curvature  of  an  isotherm  is  quite  different  for  a  two- 
phase  system.  This  difference  must  stem  from  the  variation 
in  properties  along  the  path  of  heat  flux. 

Shifting  the  temperatures,  T^  and  ,  to  moderately 
lower  values  caused  a  substantial  movement  of  the  ice/ 
water  interface.  A  useful  criterion  for  predicting  the 
interface  co-ordinates  would  be  the  ratio  of  temperature 
differences  above  and  below  the  zero  temperature  datum, 
i.e.  (T2  -  A) / (A  -  T1) ,  where  A  is  the  melting  temperature. 

4.5-3  Two-Dimensional  Isotherms 

As  previously  stated,  the  temperatures  at  discrete 
points  were  used  to  find  the  approximate  location  of  iso¬ 
thermal  lines.  The  resulting  isotherms  are  plotted  in 
fig.  4.10.  The  "geothermal"  gradient  is  noted  to  have  a 
strong  influence  on  the  isotherms  in  the  water  region. 

In  the  ice  region,  the  isotherms  tend  to  follow  the  shape 
of  the  ice/water  interface  (+32  F  isotherm) .  The  tempera¬ 
ture  gradient  across  the  insulation  panel  is  clearly  revealed. 
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Horizontal  (Inches) 

Interface  (1):  T,  =  19.5  °F,  T9  =  59.5°F  (Corresponding  Iso- 

1  o  o  therm  0.32) 

Interface  (2):  Tn  =  10.5  °F,  T9  =  36.5°F  (Corresponding  Iso¬ 
therm  0.83) 

FIG.  4.9  GRAPHICAL  DETERMINATION  OF  32 °F  ISOTHERM  VERSUS 
ACTUAL  INTERFACE  LOCATION 
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FIG.  4.10  STEADY  TWO-DIMENSIONAL  ISOTHERMS 


CHAPTER  V 


CONCLUSIONS  AND  RECOMMENDATIONS 


S.l  CONCLUSIONS 

A  simplified  theoretical  approach  to  the  problem 
of  heat  conduction  with  phase  change  has  been  presented 
in  this  thesis.  The  agreement  between  theoretical  and 
experimental  results,  for  an  ice-water  system,  using  a 
power  law  variation  of  boundary  temperature  suggests  that 
the  idealized  theory,  as  presented,  is  adequate  for  pre¬ 
dicting  ice  depths  and  temperature  profiles.  In  particu¬ 
lar,  the  assumption  of  constant  ice  properties  and  neg¬ 
lect  of  density  differences  between  the  phases  appears 
to  be  valid  for  an  ice-water  system  within  the  limitations 
of  experimental  error. 

Supercooling  of  the  liquid  phase  prior  to  formation 
of  a  solid  phase  will  occur  in  virtually  every  practical 
problem  but  has  not  been  treated  directly  in  the  ideal¬ 
ized  analytic  model.  However,  if  the  temporal  origin 
is  defined  as  that  time  when  the  temperature  is  initially 
reduced  below  the  freezing  datum,  then  the  idealized 
theory,  as  presented,  will  predict  the  solid  phase  depth 
with  reasonably  good  accuracy  once  a  stable  interface  is 
formed. 

With  the  latter  condition  realized,  crystal  growth 
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proceeds  uniformly  with  the  interface  temperature  re¬ 
maining  fixed  at  the  stable  equilibrium  freezing  point. 
If  ice  crystals  are  present  in  the  water  domain,  as 
is  evidently  the  case  for  the  periodic  surface  tempera¬ 
ture  variation  (after  the  first  cycle) ,  then  no  dis- 
cernable  supercooling  will  occur  when  the  surface  temp¬ 
erature  is  lowered  below  32 °F. 

An  approximate  power  series  solution  has  been  pre¬ 
sented  for  prediction  of  ice  depth  and  domain  tempera¬ 
ture  using  the  sinusoidual  boundary  temperature  with 
the  mean  temperature  at  the  freezing  point.  Experiment¬ 
al  results  indicate  that  the  solution  provides  a  fair 
prediction  of  the  depth  of  freezing  and  melting  for 
the  first  complete  cycle.  Theoretically  and  experiment¬ 
ally,  it  was  shown  that  a  "suspended"  slab  of  the  solid 
phase  (ice)  will  be  present  after  the  first  complete 
cycle  and  will  be  additive  to  the  overall  ice  depth  for 
the  second  complete  cycle.  If  a  positive  temperature 
gradient  exists  in  the  liquid  domain,  then,  the  maximum 
ice  depth  will  be  reduced.  The  analytic  solution,  using 
the  sinusoidual  boundary  temperature,  may  be  expected 
to  become  progressively  inaccurate  after  many  cycles, 
as  the  ice  thickness  increases.  This  occurs  because  the 
initial  domain  temperature  will  not  be  zero  at  the  be¬ 
ginning  of  each  cycle. 
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The  two-dimensional  steady  state  experiments,  using 
two  adjoining  flat  surfaces  at  different  (but  uniform) 
temperatures  above  and  below  freezing  respectively,  re¬ 
veal  that  the  interface  profile  is  quite  different  in 
curvature  from  that  predicted  graphically  for  the  cor¬ 
responding  isotherm  in  a  single  phase  system.  This 
difference  must  result  from  the  variation  of  thermal 
properties  along  the  heat  flow  path  in  a  two-phase  system. 
The  shape  of  the  interface  appears  to  be  more  sensitive 
to  any  change  in  the  temperature  of  the  surface  above 
freezing.  However,  the  co-ordinates  of  the  two-dimen¬ 
sional  steady  profile  possibly  may  be  predicted  solely 
by  the  ratio  of  surface  temperatures  above  and  below 
the  freezing  datum. 

Based  on  the  experimental  evidence  using  an  ice- 
water  system,  the  test  cell  was  larger  than  required 
for  the  particular  investigation. 


5.2  RECOMMENDATIONS 

The  simplified  theory  considers  a  domain  which  is 
initially  at  the  fusion  temperature  only.  The  dimen¬ 
sionless  formulation  can  be  extended  to  include  a  uni¬ 
form  initial  temperature  higher  than  the  fusion  tempera¬ 
ture.  The  case  of  a  constant  initial  temperature  gradient 
has  not  been  investigated.  The  use  of  computers  should 
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ultimately  aid  in  the  solution  of  the  generalized  time 
dependent  boundary  conditions  for  one-,  two-,  and  three- 
dimensional  work.  In  the  meantime,  there  is  a  pressing 
need  to  extend  the  analytic  solutions  for  practical 
boundary  conditions  as  was  done  in  the  thesis  for  a  sin- 
usoidual  variation. 

Further  experimental  tests  are  required  to  confirm 
the  simplified  theory  in  such  fields  as  casting  of  metals 
and  frost  penetration  studies.  For  laboratory  experiments 
with  solidifying  metals,  the  interface  cannot  be  physically 
observed  for  depth  measurements  as  was  the  case  for  an 
ice-water  system.  A  closed,  opaque  container  is  re¬ 
quired  for  the  study  of  solidification  of  metals;  as 
an  example,  mercury  with  its  low  melting  point  would 
require  a  relatively  simple  apparatus  and  consideration 
should  be  given  to  measuring  the  position  of  the  inter¬ 
face  on  the  basis  of  temperature  only. 
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APPENDIX  A 

FORTRAN  IV  PROGRAMME  FOR  QUADRATURE  SOLUTION 

The  Fortran  IV  programme,  for  ice  formation,  follows 
this  introduction.  Assuming  numerical  values  for  (3(.=  L) 
and  n(aD(J)),  the  programme  calculates  corresponding 
values  for  the  following,  versus  r\  =  X(I)  : 


fn  (r])  =  FU1  (I) 
o'  (ti)  =  VAR4  (I) 
Q  (t))  =  SAR3  (I) 
T  S  A8 
Ctn=-  A9 


The  particular  constant. 


’’VjrL 


*  522.0  =T,  has  been 


used  for  ice.  Changing  this  value  only  will  allow  cal¬ 


culations  for  steel,  aluminum,  and  other  materials. 
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DIMENSION  D( 15) »X( 4001 >VAR (400)  ,AR(4C0)  »SARl(400)  9  ERF ( 40  0 ) 
DIMENSION  F  U  (  4  0  o' )  9  VARl  (4  JO)  >VAR2(4D0)  >SAR2(400)  >VAR8(400)  »AKK400) 
DIMENSION  VAK4( 400 )  >SAR3 ( 400 )  >AR2 ( 400 )  > V A R 9 ( 4 0 0 )  > VAR5 ( 400 ) 
DIMENSION  VAR6 (400 ) >VAR7 (400) , DIF (400)  *FU1 (400) 

ICO  FORMAT(3X»F5«2>3X>F7«4>3X?F4»l?3X?F8*S»3X?F8*3>3X>F8*S»3X»F6«l) 

101  FORMAT  (  5X  »  1 HX  9  7 X  9  2 HF  X  9  7X  9 1 HN  9  7 X  9  3  MDU  X  9  9X  »  ^MuX  9  8  X  9  4MT  I  MIL  9  8  X  >  4nT  c.RiR  ) 
WRT  TE ( 6  » 1 0 1 ) 

L  =  4 1 
C=0.01 
C 1  =  0 . 2 
T=  52  2 . 0 
J  =  1 
D  (  JD  =  0 . 2 
K=  1 

X  (  K  )  =  0 . 0 

1  X ( K+l )  =  X ( K )  +  C 
B=~ ( X ( K ) **2 • 0 ) 

2  VAR ( K) =EXP ( B ) 

IF  (VAR(K)  .LE.  0.00001)  GO  TO  3 
K=  K+ 1 
GOT  0  1 

3  K=  1 

S AR  1  ( K )  =  0.0 

4  AR ( K )  =  (VAR(K)  +  VAR ( K+l ) )  *  C/2.0 
SARI  (K+l  )  h  SARI  CIO  +  ARCO 

5  ERF(K)  =  1.1283791  *  SARI ( K ) 

IF  (ERF(K)  .GE.  0.99998)  GO  TO  40 
K  =  K  + 1 
GOTO  4 
4  0  K=  1 

41  VAR 8 ( K )  =  ERF ( K ) /ERF ( L ) 

VA  R  9 ( K )  =  1 • U-VAR8 ( K ) 

IF  ( K  . LQ .  L )  GO  TO  6 
K  =  K  +  1 
GO  TO  41 

6  1  =  1 

X  (  I  )  =  0. 0 

7  X (  1+1 )  =  X (  I  )  +  C 
31  =  - ( X ( I )**2.5 ) 

FU (  I )  =  EXP ( B  1  ) 

8  ’  VARlt  I  )  =  EXP ( X ( I ) **2 • 0 ) 

9  VAR2 (  I  )  =  FU ( I  )  *  VARl (  I  ) 

IF  (VAR2U)  .LE.  O.oOGl)  GO  TO  10 
1  =  1  +  1 
GO  TO  7 

10  1  =  1 

ISAR2 ( I )  =  0.0 

11  AR 1(1)  =  ( VAR2 (  I  )  +  VAR2 (  1+1 )  )*C/2.0 
S A R 2  (  I  + 1  )  =  SAR2  (  I  )  +  A R IT'D' 

IF  (AR1(I)  .LE.  C. 00 001)  GO  TO  12 
1  =  1+1 
GO  TO  11 
12  1  =  1 

13  VAR4(I)  =  SAR2(1)  *  VAR(I) 

IF  (VAR(I)  .LE.  0.00001)  GO  TO  15 
1  =  1  +  1 
GO  TO  13 
15  1  =  1 

SAR3 ( I )  =  0.0 

16  A R 2 (  I )  =  ( VAR4 ( 1 )  +  VAR4 ( I +1 )  )  *  C/2.0 


. 


. 


- 


SAR3  U  +  l )  =  S A R 3  (  I  )  +  AR2  (  I  ) 

IF  ( AR2 ( I )  .LE.  0.00001)  GO  TO  34 
1  =  1  +  1 
GO  TO  16 
34  K  =  1 

33  VAR5(K)  =  VAR8(K)  *  SAR3(L) 

IF  ( K  .LQ.  L)  GO  TO  19 
K  =  K+l 
GO  TO  35 
19  1  =  1 

K  =  1 

21  VAR  6 (  I  )  =  ( SAR3 ( I )  -  VAR5(K))  *  4.0  *  D(J) 

22  FDTTrn  V A R 9  (  K  )  +  VAR6  T  I  ) 

IF  (FUKI)  .LE.  0.0001)  GO  TO  23 
B 3  =  VAR2 (  I )  /  VAR1 (  I  ) 

D I F (  I )  =  B 3  -  FU1 ( I ) 

1  =  1  +  1 
K  =  K  +  1 
GO  TO  21 

23  DO  24  1=1, L 

IF  (DIF(I)  .GT.  0.0001)  uU  TO  26 

24  CONTINUE 
A 1  =  X  (  L  ) 

A2=A1**2» 0 
A3  =  E XP ( -A 2  ) 

A4= A3 /ERF ( L) 

A5  = ( 1 . 0+4. 0*D ( J ) *SAR3 (L)  )/ ( 2 • 0*D ( J ) +1 . 0 ) 

A6=  (  0.8  86  2  26  9*4. 0*D  (  J  )*VAR4  (L)  )  /  (  2 . 0*1)  (  J  )  +1 . 0  ) 

A7  =  A 1 / ( A4*A5-A6 ) 

A  8  =  A  7  *  *  (  1  •  0  /  D  (  J  )  ) 

A9=A7*T 
I  =  1 
K=  1 

25  WR  I  T  E  (  6  >  1  0  G  )  X ( I )  ,FU1(I)  ,D(J)  , VAR 4(1 )  >  SAR3 ( I)  ,A8,A9 
1  =  1  +  1 

K  =  K+1 

IF  (I  . EQ.  ( L+lT)  GO  TO  28 
GO  TO  25 

26  1  =  1 

27  VAR7  C  I  )  =  FUKI)  *  VARKI) 

VAR2 (  I  )  =  VAR7 ( 1 ) 

IF  ( VAR2 ( I )  .LE.  0.0001 )  GO  TO  10 
I~  =1  +  1 
GO  TO  27 

28  U=  U+ 1 

D  (  J ) =D ( J- 1 ) +C 1 
I F ( D ( U  )  . GE.  3.2 )  GO  TO  36 

GO  TO  19 
36  L=L-5 

U  =  1 
D  (  U  )  =  0 . 2 

IF  ( L  . EQ.  1 )  GO  TO  30 
31  K  =  1 

33  VAR  8 ( K )  =  ERF ( K ) /ERF ( L ) 

VAR 9 ( K )  =  1.0  -  VAR8 CO 
IF  (  K  .  EQ.  L  )  -GO  TO  6 


. 
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APPENDIX  B 

COMPARISON  OF  QUADRATURE  AND  EXPLICIT  SOLUTION 


The  quadrature  and  explicit  solutions  can  be  com¬ 
pared  by  using  both  methods  in  the  numerical  evaluation 

i 

of  fn  (f3)  in  the  basic  equation  (28)  .  The  quadrature 
solution  for  the  temperature  gradient  at  £  is 


fn  (p) 


4n  q'(P)  -  -= 
Vtt 


erfp 


1  +  4n  Q(3) 


and  the  explicit  solution  is 


fn' (p) 


-22n  r(n+l)  e~ft 

Pn($)  i2n  erfc  p 


The  evaluation  shall  be  undertaken  for  integer  values 
2n  =  0,  1,  2. 


p 

1 

Quadrature  fn  (f3) 

Explicit  fn  (3) 

2n  =  0 

1 

2 

2n  = 

0 

1 

2 

o 

o 

0 

-  oo 

-  oo 

-  oo 

-  oo 

-  oo 

-  oo 

.  05 

(-) 19. 987 

(-) 19. 979 

(->19.971 

(->19. 

987 

(-) 19. 980 

(-) 19. 969 

0 

M 

O 

9.930 

9.902 

9.872 

9. 

930 

9.900 

9.868 

.15 

6.568 

6.522 

6.477 

6 . 

568 

6.519 

6.472 

o 

CM 

9 

4.868 

4.806 

4.745 

4. 

868 

4.804 

4.741 

.25 

3.836 

3.759 

3.684 

3. 

836 

3.758 

3.681 

.30 

3.138 

3.047 

2.961 

3. 

138 

3.046 

2.958 
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A  comparison  of  p  versus  fn  (|3)  shows  that 

i 

fn  (p)  ^  -  1/p  for  all  n.  Substitution  into  equation 

(28)  provides  a  second  approximation  for  dimensionless 
ice  depth,  i,e. , 


(3 


Vlr 

2 (2n+l) 


1/2 


Tn/2 


u  -  ••  c  ;  L  at.  *r  >  ijaciul  ,»  Li>  io*  \  {  i  . 
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APPENDIX  C 

LAPLACE  TRANSFORMATION  METHOD 


The  Laplace  transformation  has  been  applied  with 
considerable  success  to  problems  in  heat  conduction. 

The  purpose  of  this  section  is  to  develop  a  general  solu¬ 
tion  for  the  conduction  equation  with  reference  to  the 
one-dimensional  change  of  phase  problem.  A  solution  is 
sought  for  the  differential  equation  (1)  with  boundary 
conditions 

t  =  0  0  <  x  <  co  ©  =  ©i 

t  >  0  x  =  0  9  =  9  (t) 

T,  ° 

t  >  0  x  =  b  6  =  0 

V 

The  Laplace  transformation  of  equation  (1)  is 


|  5 


0 


and  the  general  solution  for  the  transformed  differential 
equation  is 

—  6i 

9  =  A  sinh  mx  +  B  cosh  mx  +  — 

where  S  is  the  Laplace  transform  variable  and  m  =  Ajs/ic. 
The  transformed  boundary  condition  at  x  =  0  is  9  =  0o(t). 
Substitution  of  transformed  boundary  conditions  into 
equation  (Cl)  yields 


(Cl) 


tfC  t?Av 


A  l  r  jq/>  rr^dcC  a-r  •.  i  9  1  .  n  •• 


••••:  o  * 

r  '■  «  co 1  *.  -  *:  "  3  ii  .  •.  s.-  ’ 

.-•3  >r,rrs36:- ~>:t  n3  ..->  n oj  >etp#  «o  s  r;oo  t  70  :  1  a  ... 


A 


e 
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9  =  9  (t) 


Sinh  m(b-x) 

9. 

l 

Sinh  m(b-x)  +  Sinh  mx 

S inh  mb 

S 

Sinh  mb 

+ 


(C2) 


Considering  the  initial  condition  where  9^  =  0,  equation 
(C2)  reduces  to 


9 


a  sinh  m(b-x) 

ol  '  sinh  mb 


(C3) 


The  inversion  of  equation  (C3)  is  presented  here  for  the 
simple  case  of  a  linear  variation  in  surface  temperature. 

Considering  the  temperature  variation,  9q  =  Ct, 
the  transformed  temperature  is 


9 


-mx  ,  -m(x+2b) 

e  +  e 


+  e-m(x+4b)  + 


(C4) 


Inverting  (C4)  term  by  term  and  substituting  dimensionless 
variables  where  possible,  gives 


9 


4Ct 


i^erfcT]  +  i2erfc  (r)+2f3)  +  i2erfc  (Tyt4f3)  + 


(C5 ) 


Using  the  energy  balance  equation  (16)  and  differentiating 
(C5)  term  by  term,  the  dimensionless  depth,  (3,  is  found 


to  be 

P 


^ierfcf3  +  ierfc3£  +  ierfc5(3  + 


(C6) 


The  "closed"  form  similarity  solution  for  depth  (20) , 
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using  n  ~  1,  is 


0 


-2p2  -  4i2erfcp  +  1  - 


(C7 ) 


A  comparison  of  (C6)  and  (C7)  is  made  in  the  table  fol¬ 
lowing?  consideration  need  only  be  given  to  the  numerical 
evaluation  of  the  terms  in  the  brackets.  It  can  be  seen 
that  the  series  solution  (C6)  requires  more  terms  for 
better  accuracy. 


f3 

(C6)  three  terms 
■*  «■ 

(C7)  term 

0.00 

oo 

oo 

0.05 

2.2289 

8.8485 

0.10 

1.7436 

4.3728 

0.15 

1.3410 

2.8677 

0.20 

1.0509 

2.1010 

0.25 

0.8437 

1.6312 

0. 30 

0.6930 

1.3108 

. 
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APPENDIX  D 

A  DISCUSSION  ON  DIMENSIONLESS  DEPTH,  ft  . 


As  stated  in  Chapter  II  and  shown  on  fig.  2.5,  the 
dimensionless  ice  depth  tends  to  the  limit  of  ft  =  0.765 
and  this  corresponds  roughly  to  a  surface  temperature, 
for  all  n,  of  9  (t)  =  -492 °F.  To  illustrate,  consider 
the  ice  property,  T  -  VtF  L/Cp  or 


T  o  =  522  (units  °F) 

u . 


For  the  maximum  temperature  step  change,  n  =  0,  and 
C  =  492 ,  it  is  found  that 


P 


max 


0.945 


-P‘ 


erfft 


7Z  1.0 


erfft 


8  0.765.  The  case  for  0  (t)  =  Ctn  will 

Kmax  ^  o 


or 


be  identical  for  calculating  the  limiting  value,  Pmax- 

As  a  comparison  with  other  materials,  the  following 
table  shows  the  approximate  relation  between  the 
melting  temperature  and  the  constant,  T: 


* 


Used  in  Appendix  A  for  the  computer  programme. 


£( 


,  •  *>  n  r*sctis  b  .3  X  •  -  mIO  :  i  .  n.  "  • 


S  )  ■  ■  r  •-.;■•  1  *  (  X©  ' 


1  . , a  Qi  vi  :va  i  ...bno<  id-  o  Jir  *  bn#.- 

ti  ••>  m  ■>  -  •  ‘  1  s.  /XI  : 


■ 


'  "  ""  **’ 

a  •  r  >  •  ;•  .  r  :i 

■ 

_  _ 


ifiitt  bm/ol  p.X  PX  ,  SU'fr  *  0 


i  tv  “•>  .  )  <  -rot  r.-  .f  2.T.0  ~ 


.  :■  ,  •ulcv  vJt  ti  I  ort  Jti  >fcXtpI«  >  id?  Xoa  "r  ©bi  ©cf 


£;i.  o.l  .>*  .  ^  -:  *;  irf.o  rrtf  •:  :  1  ...  t.\ 


■ 


.  ?  :  i's^aqrmo  bit)  xoi  A  >iir.  nl  bipU 


84 


Material 

Approx.  Melting 

Temp.  (  R) 

Value  of 

T  (handbook) 

Ice 

492 

522 

Aluminum 

1680 

1660 

Copper 

2440 

1910 

Iron 

3262 

3610 

Lead 

1080 

695 

Mercury 

422 

264 

Silver 

2220 

1540 

Since  a  considerable  variation  in  thermal  properties  may 
be  found  in  various  handbooks  on  materials ,  there  is  some 
reason  to  believe  that  there  should  be  actually  a  direct 
correspondence  between  melting  temperature  and  the  value 
of  T„ 

From  the  tabulation,  it  may  be  seen  that  the  constant, 
a,  as  defined  in  Chapter  II,  is  of  small  order  for  mat¬ 
erials  other  than  ice.  On  this  basis,  a  number  of  small 
order  terms  in  equation  (22)  were  neglected.  From  (22) , 
the  series  solution. 


P 


2a  tn  fn  (g) 

2n  +  1 


ft2 t1/2 

n 


B2(B2-1)  t  1 

n(n  -  1/2) 


/ 


was  developed.  The  difficulty  with  this  solution  is  that 
singularities  occur  at  n  —  0,  1/2,  1,  ...  .  However,  for 
n  -  0,  the  exact  solution  is  known  to  be  3  =  2at  f  (|3) 
and  therefore  only  the  first  term  in  the  series  can  be 
valid.  The  complete  series  solution  appears  to  be  valid 
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for  non-integer  values  of  n  or  for  long  time,  t„  Al¬ 
though  not  based  on  mathematical  rigor  but  rather  on 
the  result  from  the  case  for  n  -  0,  a  solution  seems 
possible  if  only  those  terms  prior  to  the  term  contain¬ 
ing  the  singularity  are  used.  On  this  basis,  for 
n  =  1 ,  the  first  three  terms  of  equation  may  be  used. 
Their  relative  values  are  shown  in  the  table  below, 

where  an  initial  value  for  6=6  is  assumed  and  time, 

o 

t,  calculated  using  the  first  term  of  the  series.  It 
is  noted  that  the  second  and  third  terms  are  much  smaller 
than  one  and  therefore,  the  first  term  of  the  series 
is  a  good  approximation  for  the  evaluation  of  p. 


n  =  1 

a  =  -.( 

30085  (ice) 

Assume 

Po 

p 

Q 

t 

p2  t'1/2 
Ko 

p20(p2  -  Df1 

fn  (Pn) 

n 

n(n  -  1/2) 

0 

o 

o 

.00000 

0 

.05 

.00250 

4.41 

+.00119 

-.00001 

.10 

.01013 

17.88 

+.00236 

-.00119 

.15 

.02316 

40.87 

+.00352 

-.00190 

.20 

.04215 

74.38 

+.00463 

-.00108 

.25 

.06786 

119.75 

+.00571 

-.00144 

.30 

.10132 

178.80 

+.00673 
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APPENDIX  E 

PROPERTIES  OF  TAP  WATER. 

Ordinary  tap  water,  from  the  local  region,  was  used 
in  the  experimental  work.  Typical  values  of  chemical  com¬ 
ponents  are  shown  below  for  the  month  of  January: 


Component 

River 

Treated  (Tap) 

Hardness 

241  (ppm) 

88  (ppm) 

Calcium 

195 

52 

Magnesium 

46 

36 

Alkalinity 

197 

29 

Non -Carbonates 

44 

— 

Ph  7.9  9.2 

Total  Solids  247  - 

Suspended  Solids 


24 
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